import argparse
import os.path
import subprocess
import shutil
import pandas as pd


def createDirWhenNoExist(dirPath):
    if not os.path.exists(dirPath):
        os.mkdir(dirPath)


def deleteFileIfExist(filePath):
    if os.path.exists(filePath):
        os.remove(filePath)


def deleteDirIfExist(dirPath):
    if os.path.exists(dirPath):
        shutil.rmtree(dirPath)


parser = argparse.ArgumentParser()
parser.add_argument("-l", type=str, default="", required=True, help="fq list file")
parser.add_argument("--host", type=str, default="", required=True,
                    help="* NA or species bwa index;NA mean remove human;must input")
parser.add_argument("--qual", type=int, default=33,
                    help="base quality, default=33 [Illumina is using 64, Sanger Institute is using 33](33,64)")
parser.add_argument("--contig", type=int, default=500,
                    help="min_contig_length,default: 500")
parser.add_argument("--identify", type=int, default=0.95,
                    help="cd-hit for gene,identify,default: 0.95")
parser.add_argument("--overlap", type=int, default=0.9,
                    help="cd-hit for gene,overlap,default: 0.9")
args = parser.parse_args()

scriptDir = os.path.dirname(os.path.abspath(__file__))
listFile = args.l
qual = args.qual
contig = args.contig
identify = args.identify
overlap = args.overlap
ref = "00.ref"
createDirWhenNoExist(ref)
host = args.host
binDir = scriptDir + "/bin"
databaseDir = scriptDir + "/database"
pwd = os.getcwd()
qc_fi = "01.QC"
gene = "05.geneset"
createDirWhenNoExist(gene)

# generate qc script
cmd = f"perl {binDir}/QC_work.pl {listFile} -qual {qual}"
print(cmd)
subprocess.check_call(cmd, shell=True)

# generate index script
with open(f"shell/S00.ref.index.sh", "w") as file:
    file.write(f"cd {ref}/\n")
    if host == "NA":
        file.write(f"ln -s {databaseDir}/host_db/human.fasta.amb host.fa.amb\n")
        file.write(f"ln -s {databaseDir}/host_db/human.fasta.ann host.fa.ann\n")
        file.write(f"ln -s {databaseDir}/host_db/human.fasta.bwt host.fa.bwt\n")
        file.write(f"ln -s {databaseDir}/host_db/human.fasta.pac host.fa.pac\n")
        file.write(f"ln -s {databaseDir}/host_db/human.fasta.sa  host.fa.sa\n")
    else:
        file.write(f"cat {binDir}/host_db/human.fasta {host} >host.fa\n")
        file.write(f"{binDir}/software/bwa index host.fa\n")
    file.close()

# produce S01 script
asse = "03.Assembly"
createDirWhenNoExist(asse)
asseVir = "03.Assembly.vir"
createDirWhenNoExist(asseVir)
loca = "03.Assembly.loca"
createDirWhenNoExist(loca)
pred = "04.predict"
createDirWhenNoExist(pred)
sampleNum = sum(1 for _ in open(listFile))
locaLstFile = f"{pwd}/{loca}/loca.lst"
predictGeneListFile = f"{pwd}/predict.gene.list"
hostInfoLstFile = f"{pwd}/{qc_fi}/host.info.lst"
prof = "06.profile"
createDirWhenNoExist(prof)
samListFile = f"{pwd}/{prof}/sam.lst"
viralSamListFile = f"{pwd}/{prof}/viral.sam.lst"
shell = "shell"
createDirWhenNoExist(shell)
with open(f"{listFile}", "r") as file:
    locaLstFileH = open(locaLstFile, "w")
    predictGeneListFileH = open(predictGeneListFile, "w")
    hostInfoLstFileH = open(hostInfoLstFile, "w")
    samListFileH = open(samListFile, "w")
    viralSamListFileH = open(viralSamListFile, "w")
    for line in file:
        columns = line.strip().split("\t")
        sampleName = columns[0]
        outFileName = f"shell/S01.1.{sampleName}.sh"
        outFile = open(outFileName, "r")
        outFile.write(
            f"{binDir}/software/bwa mem -t 40 -M {pwd}/00.ref/host.fa {sampleName}.clip.1.fq.gz {sampleName}.clip.2.fq.gz |awk '$3!=\"*\"'  > {sampleName}.host.sam\n")
        outFile.write(f"perl {binDir}/Meta/remove-host.pl {sampleName}.host.sam {sampleName}\n")
        outFile.write(f"rm {sampleName}.host.sam\n")
        hostInfoLstFileH.write(f"{sampleName}\t{pwd}/{qc_fi}/{sampleName}.host.info\n")
        outFile.write(f"\ncd {pwd}/{asse}/\n")
        outFile.write(
            f"{binDir}/software/megahit/megahit -1 {pwd}/{qc_fi}/{sampleName}.clean.1.fq.gz -2 {pwd}/{qc_fi}/{sampleName}.clean.2.fq.gz --min-contig-len {contig} -t 40 -o {pwd}/{asse}/{sampleName}\n")
        outFile.write(
            f"perl {binDir}/Meta/renamefa.pl {pwd}/{asse}/{sampleName}/final.contigs.fa {sampleName} {pwd}/{asse}/{sampleName}.contig.ok.fa\n")
        outFile.write(
            f"perl {binDir}/deal_fa.pl -format 3 {pwd}/{asse}/{sampleName}.contig.ok.fa |perl -e  'while(<>){{chomp;@a=split; if($a[1] > 10000){{$a[1]=10000;}} print \"$a[0]\\t$a[1]\\n\";}}' >{pwd}/{asse}/{sampleName}.contig.ok.fa.chrlist\n")
        outFile.write(
            f"perl {binDir}/fa_fq_len_bar.pl {pwd}/{asse}/{sampleName}.contig.ok.fa.chrlist {pwd}/{asse}/{sampleName}.contig.length.pdf contig\n")
        outFile.write(f"rm -r {pwd}/{asse}/{sampleName}/\n\n")
        outFile.write(f"cd {pwd}/{loca}/\n")
        outFile.write(
            f"perl {binDir}/deal_fa.pl {pwd}/{asse}/{sampleName}.contig.ok.fa -format 6 -len 1000 >{sampleName}.contig.1k\n")
        outFile.write(
            f"perl {binDir}/deal_fa.pl {pwd}/{asse}/{sampleName}.contig.ok.fa -len 5000 -format 6 -type 1 | perl {binDir}/deal_fa.pl - -len 1000000 -format 6 -type 2 >{sampleName}.contig.ok.5k\n")
        outFile.write(f"perl {binDir}/deal_fa.pl {sampleName}.contig.1k {sampleName}.split -format 11 -len 40000\n")
        outFile.write(f"source activate plasflow\n")
        outFile.write(f"ls {sampleName}.split*fa |while read dd;do PlasFlow.py --input $dd --output $dd.out;done\n")
        outFile.write(f"conda deactivate\n")
        outFile.write(f"cat {sampleName}.split*fa.out |cut -f 3-         >{sampleName}.location.xls\n")
        outFile.write(f"cat {sampleName}.split*fa.out_chromosomes.fasta  >{sampleName}.location.chromosome.fa\n")
        outFile.write(f"cat {sampleName}.split*fa.out_plasmids.fasta     >{sampleName}.location.plasmid.fa\n")
        outFile.write(f"cat {sampleName}.split*fa.out_unclassified.fasta >{sampleName}.location.unclassified.fa\n")
        outFile.write(f"rm {sampleName}.split*fa* {sampleName}.contig.1k\n\n")
        locaLstFileH.write(
            f"{sampleName}\t{pwd}/{loca}/{sampleName}.location.chromosome.fa\t{pwd}/{loca}/{sampleName}.location.plasmid.fa\t{pwd}/{loca}/{sampleName}.location.unclassified.fa\n")
        outFile.write(f"cd {pwd}/{pred}/\n")
        outFile.write(
            f"{binDir}/software/prodigal -a {sampleName}.contig.orf.faa -i {pwd}/{asse}/{sampleName}.contig.ok.fa -o {sampleName}.temp.txt -p meta -q -d {sampleName}.contig.orf.ffn\n")
        if sampleNum > 6:
            outFile.write(
                f"awk '$2 >= 1500' {pwd}/{asse}/{sampleName}.contig.ok.fa.chrlist  |perl {binDir}/Meta/extract.contig.gene.pl - {sampleName}.contig.orf.ffn >{sampleName}.contig.orf.ffn.fil\n")
        else:
            outFile.write(f"ln -s {sampleName}.contig.orf.ffn {sampleName}.contig.orf.ffn.fil\n")
        outFile.write(f"gzip -f {pwd}/{asse}/{sampleName}.contig.ok.fa\n")
        outFile.write(f"rm {sampleName}.temp.txt\n")
        predictGeneListFileH.write(f"{sampleName}\t{pwd}/{pred}/{sampleName}.contig.orf.ffn\n")
        outFile.write(
            f"{binDir}/software/cd-hit/cd-hit-est -i {sampleName}.contig.orf.ffn.fil -o {pwd}/{gene}/{sampleName}.geneSet.ffn -n 9 -c {identify} -G 0 -M 0 -d 0 -aS {overlap} -r 0 -T 40\n")
        outFile.write(f"rm {sampleName}.contig.orf.ffn.fil\n")
        outFile.close()
        mapFile = f"{pwd}/{shell}/S02.1.{sampleName}.map.sh"
        mapFileH = open(mapFile, "w")
        mapFileH.write(f"gzip -f {pwd}/{pred}/{sampleName}.contig.orf.faa\n")
        mapFileH.write(f"gzip -f {pwd}/{pred}/{sampleName}.contig.orf.ffn\n")
        mapFileH.write(f"cd {pwd}/{gene}/\n")
        mapFileH.write(
            f"{binDir}/software/bwa mem -t 40 -M -R '@RG\tID:{sampleName}\tSM:{sampleName}\tLB:{sampleName}\tPL:Illumina\tPI:500' uniqGeneSet.fa {pwd}/{qc_fi}/{sampleName}.clean.1.fq.gz {pwd}/{qc_fi}/{sampleName}.clean.2.fq.gz > {sampleName}.sam\n")
        mapFileH.write(
            f"perl {binDir}/Meta/uniqGene-ProfilebyBwa.pl uniqGeneSet.fa {sampleName}.sam {sampleName}.sam.abd\n")
        samListFileH.write(f"{sampleName}\t{pwd}/{gene}/{sampleName}.sam.abd\n")
        mapFileH.write(f"cd {pwd}/{prof}/\n")
        mapFileH.write(
            f"{binDir}/software/bwa mem -t 40 -M -R '@RG\\tID:{sampleName}\\tSM:{sampleName}\\tLB:{sampleName}\\tPL:Illumina\\tPI:500' {pwd}/{pred}/phagcn2.viral.fa {pwd}/{qc_fi}/{sampleName}.clean.1.fq.gz {pwd}/{qc_fi}/{sampleName}.clean.2.fq.gz > {sampleName}.viral.sam\n")
        mapFileH.write(
            f"perl {binDir}/Meta/uniqGene-ProfilebyBwa.pl {pwd}/{pred}/phagcn2.viral.fa {sampleName}.viral.sam {sampleName}.viral.sam.abd\n")
        mapFileH.write(f"rm {sampleName}.viral.bam {sampleName}.viral.sam\n")
        viralSamListFileH.write(f"{sampleName}\t{pwd}/{prof}/{sampleName}.viral.sam.abd\n")
        mapFileH.close()
    locaLstFileH.close()
    samListFileH.close()
    viralSamListFileH.close()
    predictGeneListFileH.close()
    hostInfoLstFileH.close()

comb = "08.ComponentAnaly"
createDirWhenNoExist(comb)
comb1 = "08.ComponentAnaly/Heatmap"
createDirWhenNoExist(comb1)
comb2 = "08.ComponentAnaly/Venn"
createDirWhenNoExist(comb2)
s01StatFile = f"{pwd}/{shell}/S01.2.stat.sh"
deleteFileIfExist(s01StatFile)
vira = "04.Viral"
createDirWhenNoExist(vira)
with open(s01StatFile, "w") as file:
    file.write(f"perl {binDir}/Meta/2-2.assembly_stat.pl   {pwd}/{asse}/ {contig} {pwd}/{asse}/\n")
    file.write(f"perl {binDir}/Meta/location.stat.pl {pwd}/{loca}/loca.lst > {pwd}/{loca}/location.stat.xls\n")
    file.write(f"perl {binDir}/Meta/stat_gene_predict.pl {pwd}/predict.gene.list > {pwd}/{pred}/predictGene.stat.xls\n")
    file.write(f"rm {pwd}/{asse}/*ok.fa.gz.stat {pwd}/{asse}/ass_contig.list {pwd}/predict.gene.list\n")
    file.write(f"\ncd {pwd}/{gene}/\n")
    file.write(f"cat *.geneSet.ffn > all.sam.geneSet.fa\n")
    file.write(
        f"{binDir}/software/cd-hit/cd-hit-est -i all.sam.geneSet.fa -o uniqGeneSet.fa -n 9 -c {identify} -G 0 -M 0 -d 0 -aS {overlap} -r 0 -T 96\n\n")
    file.write(f"{binDir}/software/bwa index uniqGeneSet.fa\n")
    file.write(f"rm {pwd}/{gene}/*.geneSet.ffn*\n\n")
    file.write(
        f"perl {binDir}/Meta/merge_host_info.pl {pwd}/{qc_fi}/host.info.lst >{pwd}/{qc_fi}/sample.host.stat.xls\n\n")
    file.write(f"\ncd {pwd}/{vira}/\n")
    file.write(f"cat {pwd}/{loca}/*ok.5k > ass.5k.fa\n")
    file.write(f"source activate deepvirfinder\n")
    file.write(f"condaPath=$(conda info --base)\n")
    file.write(
        f"python $condaPath/envs/deepvirfinder/bin/dvf.py -i {pwd}/{vira}/ass.5k.fa -o ./ -l 300 -c 64 -m $condaPath/envs/deepvirfinder/models\n")
    file.write(f"conda deactivate\n")
    file.write(
        f"awk '$3 >= 0.7 && $4 <=0.05' {pwd}/{vira}/ass.5k.fa_gt300bp_dvfpred.txt |awk 'BEGIN{{print \"ID\\tlength\\tscorel\\tpvalue\";}}{{print $0}}' > VirFinder.viral.xls\n")
    file.write(f"perl {binDir}/find_identical.pl ass.5k.fa VirFinder.viral.xls -format 2 >VirFinder.viral.fa\n\n")
    file.write(f"source activate vs2\n")
    file.write(
        f"virsorter run -w virsorter2-out --include-groups dsDNAphage,ssDNA,NCLDV -i VirFinder.viral.fa -j 96 all\n")
    file.write(f"conda deactivate\n")
    file.write(
        f"awk '$1!~/lt2gene/ && $1!~/partial/' virsorter2-out/final-viral-score.tsv |perl -e '<>; print \"ContigID\\ttype\\n\"; while(<>){{chomp;@a=split /\\t/; next if($a[4] < 0.95); @b=split /\\|\\|/,$a[0];print \"$b[0]\\t$a[5]\\n\";}}'  >virsorter.viral.xls\n")
    file.write(f"perl {binDir}/find_identical.pl ass.5k.fa virsorter.viral.xls -format 2 >virsorter.viral.fa\n\n")
    file.write(f"cd {pwd}/{pred}\n")
    file.write(f"source activate checkv\n")
    file.write(f"checkv end_to_end {pwd}/{vira}/virsorter.viral.fa ./ -t 64 -d $condaPath/envs/checkv/checkv-db-v1.4\n")
    file.write(f"conda deactivate\n\n")
    file.write(f"cut -f 1-3,6-12 quality_summary.tsv |sed -e 's/contig_id/vOTU/' >checkv.virus.summary.xls\n")
    file.write(
        f"cut -f 6 checkv.virus.summary.xls |sed -n '2,$'p |sort |uniq -c |awk 'BEGIN{{print \"type\\tnum\\tnum2\";}}{{print $2\"\\t\"$1\"\\t\"$1}}' >checkv.virus.stat.xls\n")
    file.write(f"perl {binDir}/bar_pie.pl -i checkv.virus.stat.xls -bar F -p 0\n")
    file.write(f"rm pie.num2.checkv.virus.stat.xls.pdf; mv pie.num.checkv.virus.stat.xls.pdf pie.checkv.virus.pdf\n")
    file.write(f"python3 {scriptDir}/checkv_region.py")
    file.write(f"samtools faidx {pwd}/{vira}/virsorter.viral.fa -r checkv_region.xls -o checkv.viral.fa\n")
    file.write(f"source activate bbmap\n")
    file.write(
        f"dedupe.sh in=checkv.viral.fa out=uniq.viral.fna minidentity=95 minlengthpercent=85 minoverlappercent=85 overwrite=t pbr=t cluster=t\n")
    file.write(f"conda deactivate\n\n")
    file.write(f"source activate phagcn2\n")
    file.write(f"export MKL_SERVICE_FORCE_INTEL=1\n")
    file.write(f"cd  $condaPath/envs/PhaGCN2.0-PhaGCNv2.0/\n")
    file.write(f"python run_Speed_up.py --contigs {pwd}/{pred}/uniq.viral.fna --len 8000\n")
    file.write(f"conda deactivate\n\n")
    file.write(f"cp final_prediction.csv {pwd}/{pred}/\n")
    file.write(f"cd {pwd}/{pred}\n")
    file.write(f"samtools faidx uniq.viral.fna -r phagcn2_region.xls -o phagcn2.viral.fa\n")
    file.write(f"{binDir}/software/bwa index         phagcn2.viral.fa\n")
    file.write(f"perl {binDir}/deal_fa.pl -format 3  phagcn2.viral.fa > phagcn2.viral.fa.chrlist\n\n")

s02StatFile = f"{pwd}/{shell}/S02.2.stat.sh"
deleteFileIfExist(s02StatFile)
with open(s02StatFile, "w") as file:
    file.write(f"cd {pwd}/{prof}/\n")
    file.write(f"perl {binDir}/Meta/Profile-Merge-BigBAM.pl sam.lst temp.profile.txt temp.reads.txt\n")
    file.write(
        f"less temp.profile.txt  | perl -e 'my $head = <>;print $head;while(<>){{chomp;my @a = split;my $b = 0;for ($i = 1;$i<=$#a;$i++){{$b+=$a[$i];}}if($b>0){{print \"$_\n\";}}}}' > profile.txt\n")
    file.write(f"perl {binDir}/find_identical.pl temp.reads.txt profile.txt > reads.txt\n")
    file.write(f"rm -rf sam.lst temp.profile.txt  temp.reads.txt {pwd}/{gene}/*.sam*\n")
    file.write(
        f"perl {binDir}/find_identical.pl {pwd}/{gene}/uniqGeneSet.fa profile.txt -format 2 |awk '{{print $1}}' > uniqGeneSet.ffn\n")
    file.write(f"{binDir}/transeq -sequence uniqGeneSet.ffn -table 11 -trim -outseq uniqGeneSet.faa\n")
    file.write(f"sed -i 's/_1$//' uniqGeneSet.faa\n")
    file.write(
        f"perl {binDir}/deal_fa.pl -format 3 uniqGeneSet.ffn |perl -e 'while(<>){{chomp;@a=split; if($a[1] > 5000){{$a[1]=5000;}} print \"$a[0]\\t$a[1]\\n\";}}' >uniqGeneSet.ffn.chrlist\n")
    file.write(f"perl {binDir}/fa_fq_len_bar.pl uniqGeneSet.ffn.chrlist uniqGeneSet.length.dist.pdf gene\n")
    file.write(f"rm {pwd}/{gene}/uniqGeneSet.fa* uniqGeneSet.ffn.chrlist uniqGeneSet.length.dist.pdf.cmd.r\n")
    file.write(f"cd {pwd}/{comb2}\n")
    file.write(f"perl {binDir}/exp.mat.venn.pl {pwd}/{prof}/profile.txt\n")
    file.write(f"rm *.imp cmd.r\n\n")
    file.write(f"perl {binDir}/Meta/Profile-Merge-BigBAM_virus.pl viral.sam.lst viral.profile.txt viral.reads.txt\n")

s03NrFile = f"{pwd}/{shell}/S03.nr.sh"
deleteFileIfExist(s03NrFile)
anno = "07.Annotation"
createDirWhenNoExist(anno)
nr = f"{anno}/NR"
createDirWhenNoExist(nr)
sarg = "07.Annotation/SARG"
createDirWhenNoExist(sarg)
with open(s03NrFile, "w") as file:
    file.write(f"cd {pwd}/{anno}/\n")
    file.write(
        f"{binDir}/software/diamond blastp -d {databaseDir}/gtdb.anno/00.index/gtdbtk.bar.arc -q {pwd}/{prof}/uniqGeneSet.faa -o uniqGeneSet.m8 -f 6 --evalue 0.00001 -k 10 -t ./ -b 8\n")
    file.write(f"perl {binDir}/best_m8.pl uniqGeneSet.m8 >uniqGeneSet.m8.fil\n\n")
    file.write(f"cd {pwd}/{nr}/\n")
    file.write(
        f"perl {binDir}/gtdb.anno/gtdb.anno2nr.pl   ../uniqGeneSet.m8.fil {databaseDir}/gtdb.anno/03.NR/GTDB-NR.bac.arc.Info >uniqGeneSet.NR.anno.xls")
    file.write(
        f"perl {binDir}/Meta/4-6.tax_profile.pl -i uniqGeneSet.NR.anno.xls -p {pwd}/{prof}/profile.txt -o uniqGeneSet.nr.tax.profile.xls\n")
    file.write(f"perl {binDir}/Meta/4-7.taxLevel_profile.pl -i uniqGeneSet.nr.tax.profile.xls -o summary\n")
    file.write(f"cd {pwd}/{nr}/summary\n")
    file.write(f"perl {binDir}/zhaozl-figure/meta.acc.bar.pl phylum.xls phylum.bar.pdf Phylum\n")
    file.write(f"perl {binDir}/zhaozl-figure/meta.acc.bar.pl genus.xls  genus.bar.pdf  Genus\n")
    file.write(f"source activate R422\n")
    file.write(f"perl {binDir}/zhaozl-figure/meta.anno.pie.pl phylum.xls phylum.pie.pdf\n")
    file.write(f"perl {binDir}/zhaozl-figure/meta.anno.pie.pl genus.xls  genus.pie.pdf\n")
    file.write(f"cd {pwd}/{comb1}/\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.type.pl {pwd}/{nr}/summary/phylum.xls heatmap_phylum.pdf 1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.type.pl {pwd}/{nr}/summary/genus.xls  heatmap_genus.pdf  1\n\n")
    file.write(f"conda deactivate\n\n")

s03SargFile = f"{pwd}/{shell}/S03.sarg.sh"
deleteFileIfExist(s03SargFile)
bacRegionFile = f"{pwd}/{nr}/bac.region.xls"
arcRegionFile = f"{pwd}/{nr}/arc.region.xls"
with open(s03SargFile, "w") as file:
    file.write(f"cd {pwd}/{prof}\n")
    file.write(f"samtools faidx uniqGeneSet.faa -r {bacRegionFile} -o uniqGeneSet.bac.faa\n")
    file.write(f"samtools faidx uniqGeneSet.faa -r {arcRegionFile} -o uniqGeneSet.arc.faa\n")
    file.write(f"cd {pwd}/{sarg}\n")
    file.write(
        f"{binDir}/software/diamond blastp -d {databaseDir}/Anno/SARG/SARG -q {pwd}/{prof}/uniqGeneSet.bac.faa -o sarg.bac.m8 -f 6 --evalue 0.00001 --id 70 -k 10 -t ./ -b 8\n")
    file.write(
        f"{binDir}/software/diamond blastp -d {databaseDir}/Anno/SARG/SARG -q {pwd}/{prof}/uniqGeneSet.arc.faa -o sarg.arc.m8 -f 6 --evalue 0.00001 --id 70 -k 10 -t ./ -b 8\n")
    file.write(f"cp {pwd}/{pred}/phagcn2.viral.fa {pwd}/{prof}/uniqGeneSet.vir.faa\n")
    file.write(
        f"{binDir}/software/diamond blastp -d {databaseDir}/Anno/SARG/SARG -q {pwd}/{prof}/uniqGeneSet.vir.faa -o sarg.vir.m8 -f 6 --evalue 0.00001 --id 70 -k 10 -t ./ -b 8\n")
    file.write(f"perl {binDir}/best_m8.pl sarg.bac.m8 >sarg.bac.m8.fil\n")
    file.write(f"perl {binDir}/best_m8.pl sarg.vir.m8 >sarg.vir.m8.fil\n")
    file.write(
        f"perl {binDir}/Anno/add.m8.phi.func.pl sarg.bac.m8.fil {databaseDir}/Anno/SARG/SARG.func.xls >uniqGeneSet.SARG.bac.anno.xls\n")
    file.write(
        f"perl {binDir}/Anno/add.m8.phi.func.pl sarg.arc.m8.fil {databaseDir}/Anno/SARG/SARG.func.xls >uniqGeneSet.SARG.arc.anno.xls\n")
    file.write(
        f"perl {binDir}/Anno/add.m8.phi.func.pl sarg.vir.m8.fil {databaseDir}/Anno/SARG/SARG.func.xls >uniqGeneSet.SARG.vir.anno.xls\n")
    file.write(
        f"cut -f 1,5- uniqGeneSet.SARG.bac.anno.xls |perl {binDir}/Meta/add.profile.pl {pwd}/{prof}/profile.txt - >uniqGeneSet.SARG.bac.profile.xls\n")
    file.write(
        f"cut -f 1,5- uniqGeneSet.SARG.arc.anno.xls |perl {binDir}/Meta/add.profile.pl {pwd}/{prof}/profile.txt - >uniqGeneSet.SARG.arc.profile.xls\n")
    file.write(
        f"cut -f 1,5- uniqGeneSet.SARG.vir.anno.xls |perl {binDir}/Meta/add.profile.pl {pwd}/{prof}/viral.profile.txt - >uniqGeneSet.SARG.vir.profile.xls\n")
    file.write(f"perl {binDir}/Meta/st.sarg.profile.pl uniqGeneSet.SARG.bac.profile.xls uniqGeneSet.SARG.bac\n\n")
    file.write(f"perl {binDir}/Meta/st.sarg.profile.pl uniqGeneSet.SARG.arc.profile.xls uniqGeneSet.SARG.arc\n\n")
    file.write(f"perl {binDir}/Meta/st.sarg.profile.pl uniqGeneSet.SARG.vir.profile.xls uniqGeneSet.SARG.vir\n\n")
    file.write(f"source activate R422\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.pie.pl          uniqGeneSet.SARG.bac.type.profile.xls      SARG.bac.type.profile.pie.pdf\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.pie.pl          uniqGeneSet.SARG.arc.type.profile.xls      SARG.arc.type.profile.pie.pdf\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.pie.pl          uniqGeneSet.SARG.vir.type.profile.xls      SARG.vir.type.profile.pie.pdf\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.pie.pl          uniqGeneSet.SARG.bac.Mechanism.profile.xls SARG.Mechanism.bac.profile.pie.pdf\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.pie.pl          uniqGeneSet.SARG.arc.Mechanism.profile.xls SARG.Mechanism.arc.profile.pie.pdf\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.pie.pl          uniqGeneSet.SARG.vir.Mechanism.profile.xls SARG.Mechanism.vir.profile.pie.pdf\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.pl      uniqGeneSet.SARG.bac.subtype.profile.xls   SARG.bac.subtype.profile.heatmap.pdf   1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.pl      uniqGeneSet.SARG.arc.subtype.profile.xls   SARG.arc.subtype.profile.heatmap.pdf   1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.pl      uniqGeneSet.SARG.vir.subtype.profile.xls   SARG.vir.subtype.profile.heatmap.pdf   1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.type.pl uniqGeneSet.SARG.bac.type.profile.xls      SARG.bac.type.profile.heatmap.pdf      1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.type.pl uniqGeneSet.SARG.arc.type.profile.xls      SARG.arc.type.profile.heatmap.pdf      1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.type.pl uniqGeneSet.SARG.vir.type.profile.xls      SARG.vir.type.profile.heatmap.pdf      1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.type.pl uniqGeneSet.SARG.bac.Mechanism.profile.xls SARG.bac.Mechanism.profile.heatmap.pdf 1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.type.pl uniqGeneSet.SARG.arc.Mechanism.profile.xls SARG.arc.Mechanism.profile.heatmap.pdf 1\n")
    file.write(
        f"perl {binDir}/zhaozl-figure/meta.anno.heatmap.type.pl uniqGeneSet.SARG.vir.Mechanism.profile.xls SARG.vir.Mechanism.profile.heatmap.pdf 1\n")
    file.write(f"conda deactivate\n")

files = os.listdir(ref)
print(files)


def deleteFile(x):
    os.unlink(x)


list(map(lambda x: deleteFile(f"{ref}/{x}"), files))
cmd = "sh shell/S00.ref.index.sh"
print(cmd)
subprocess.check_call(cmd, shell=True)

# run S01 script
files = os.listdir("shell")
step1ShFiles = list(filter((lambda x: x.startswith("S01.1")), files))
cmd = " && ".join(list(map((lambda x: f"sh shell/{x}"), step1ShFiles)))
print(cmd)
subprocess.check_call(cmd, shell=True)

# run S01 stat script
cmd = f"sh {s01StatFile}"
print(cmd)
subprocess.check_call(cmd, shell=True)

# run S02 script
files = os.listdir("shell")
stepShFiles = list(filter((lambda x: x.startswith("S02.1")), files))
cmd = " && ".join(list(map((lambda x: f"sh shell/{x}"), stepShFiles)))
print(cmd)
subprocess.check_call(cmd, shell=True)

# run S02 Stat script
cmd = f"sh {s02StatFile}"
print(cmd)
subprocess.check_call(cmd, shell=True)

# run S03 NR script
cmd = f"sh {s03NrFile}"
print(cmd)
subprocess.check_call(cmd, shell=True)

# generate bac arc region file
nrAnnoFile = f"{nr}/uniqGeneSet.NR.anno.xls"
df = pd.read_csv(nrAnnoFile, sep="\t")
bacDf = df[df["Tax"].str.startswith("d__Bacteria;")]["GeneID"]
arcDf = df[df["Tax"].str.startswith("d__Archaea;")]["GeneID"]
bacDf.to_csv(bacRegionFile, sep="\t", header=False, index=False)
arcDf.to_csv(arcRegionFile, sep="\t", header=False, index=False)

# run S03 SARG script
cmd = f"sh {s03SargFile}"
print(cmd)
subprocess.check_call(cmd, shell=True)
